如何使用GDAL重采样图像

最新推荐文章于 2024-04-19 03:48:37 发布
最新推荐文章于 2024-04-19 03:48:37 发布
阅读量2.5w 收藏 37
点赞数 4
分类专栏： GDAL C++编程技术 GDAL交流 文章标签： null 算法 algorithm float windows string
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/liminlu0314/article/details/6130064
版权
GDAL 同时被 3 个专栏收录
122 篇文章 144 订阅
订阅专栏
GDAL交流
114 篇文章 241 订阅
订阅专栏
C++编程技术
61 篇文章 8 订阅
订阅专栏

    在编写重采样图像时，可以使用GDAL来读写图像，然后自己编写重采样算法（最邻近像元法，双线性内插法，三次立方卷积法等）【关于这采样算法有时间我会单独写一篇文章来说明原理的】将计算的结果写入图像中来实现。

    在GDAL的算法中，已经提供了五种重采样算法，其定义如下（位置gdalwarper.h 的46行）：

/*! Warp Resampling Algorithm */
typedef enum {
  /*! Nearest neighbour (select on one input pixel) */ GRA_NearestNeighbour=0,
  /*! Bilinear (2x2 kernel) */                         GRA_Bilinear=1,
  /*! Cubic Convolution Approximation (4x4 kernel) */  GRA_Cubic=2,
  /*! Cubic B-Spline Approximation (4x4 kernel) */     GRA_CubicSpline=3,
  /*! Lanczos windowed sinc interpolation (6x6 kernel) */ GRA_Lanczos=4
} GDALResampleAlg;

    在查看Gdalwarp的源代码发现，warp的功能非常强大，可以用来做投影转换，重投影，投影定义，重采样，镶嵌，几何精校正和影像配准等。一句话，很好很强大。下面就看看其中的一点点皮毛，使用warp来编写一个重采样的接口，代码如下：

/**
* 重采样函数(GDAL)
* @param pszSrcFile        输入文件的路径
* @param pszOutFile        写入的结果图像的路径
* @param fResX             X转换采样比，默认大小为1.0，大于1图像变大，小于1表示图像缩小
* @param fResY             Y转换采样比，默认大小为1.0
* @param nResampleMode     采样模式，有五种，具体参见GDALResampleAlg定义，默认为双线性内插
* @param pExtent           采样范围，为NULL表示计算全图
* @param pBandIndex        指定的采样波段序号，为NULL表示采样全部波段
* @param pBandCount        采样的波段个数，同pBandIndex一同使用，表示采样波段的个数
* @param pszFormat         写入的结果图像的格式
* @param pProgress         进度条指针
* @return 成功返回0，否则为其他值
*/
int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, float fResX , float fResY, LT_ResampleMode nResampleMode,
    LT_Envelope* pExtent, int* pBandIndex, int *pBandCount, const char* pszFormat,  LT_Progress *pProgress)
{
    if(pProgress != NULL)
    {
        pProgress->SetProgressCaption("重?采?样?");
        pProgress->SetProgressTip("正?在?执?行?重?采?样?...");
    }

    GDALAllRegister();
    GDALDataset *pDSrc = (GDALDataset *)GDALOpen(pszSrcFile, GA_ReadOnly);
    if (pDSrc == NULL)
    {
        if(pProgress != NULL)
            pProgress->SetProgressTip("指?定?的?文?件?不?存?在?，?或?者?该?格?式?不?被?支?持?！?");

        return RE_NOFILE;
    }

    GDALDriver *pDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if (pDriver == NULL)
    {
        if(pProgress != NULL)
            pProgress->SetProgressTip("不?能?创?建?该?格?式?的?文?件?！?");

        GDALClose((GDALDatasetH) pDSrc);
        return RE_CREATEFILE;
    }

    int iBandCount = pDSrc->GetRasterCount();
    string strWkt = pDSrc->GetProjectionRef();
    GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();

    double dGeoTrans[6] = {0};
    pDSrc->GetGeoTransform(dGeoTrans);

    int iNewBandCount = iBandCount;
    if (pBandIndex != NULL && pBandCount != NULL)
    {
        int iMaxBandIndex = pBandIndex[0];    //找?出?最?大?的?波?段?索?引?序?号?
        for (int i=1; i<*pBandCount; i++)
        {
            if (iMaxBandIndex < pBandIndex[i])
                iMaxBandIndex = pBandIndex[i];
        }

        if(iMaxBandIndex > iBandCount)
        {
            if(pProgress != NULL)
                pProgress->SetProgressTip("指?定?的?波?段?序?号?超?过?图?像?的?波?段?数?，?请?检?查?输?入?参?数?！?");

            GDALClose((GDALDatasetH) pDSrc);
            return RE_PARAMERROR;
        }
        
        iNewBandCount = *pBandCount;
    }

    LT_Envelope enExtent;
    enExtent.setToNull();

    if (pExtent == NULL)    //全?图?计?算?
    {
        double dPrj[4] = {0};    //x1,x2,y1,y2
        ImageRowCol2Projection(dGeoTrans, 0, 0, dPrj[0], dPrj[2]);
        ImageRowCol2Projection(dGeoTrans, pDSrc->GetRasterXSize(), pDSrc->GetRasterYSize(), dPrj[1], dPrj[3]);
        enExtent.init(dPrj[0], dPrj[1], dPrj[2], dPrj[3]);

        pExtent = &enExtent;
    }
    
    dGeoTrans[0] = pExtent->getMinX();
    dGeoTrans[3] = pExtent->getMaxY();
    dGeoTrans[1] = dGeoTrans[1] / fResX;
    dGeoTrans[5] = dGeoTrans[5] / fResY;

    int iNewWidth  = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[1]) + 0.5) );
    int iNewHeight = static_cast<int>( (pExtent->getMaxX() - pExtent->getMinX() / ABS(dGeoTrans[5]) + 0.5) );

    GDALDataset *pDDst = pDriver->Create(pszOutFile, iNewWidth, iNewHeight, iNewBandCount, dataType, NULL);
    if (pDDst == NULL)
    {
        if(pProgress != NULL)
            pProgress->SetProgressTip("创?建?输?出?文?件?失?败?！?");

        GDALClose((GDALDatasetH) pDSrc);
        return RE_CREATEFILE;
    }

    pDDst->SetProjection(strWkt.c_str());
    pDDst->SetGeoTransform(dGeoTrans);

    GDALResampleAlg eResample = (GDALResampleAlg) nResampleMode;

    if(pProgress != NULL)
    {
        pProgress->SetProgressTip("正?在?执?行?重?采?样?...");
        pProgress->SetProgressTotalStep(iNewBandCount*iNewHeight);
    }

    int *pSrcBand = NULL;
    int *pDstBand = NULL;
    int iBandSize = 0;
    if (pBandIndex != NULL && pBandCount != NULL)
    {
        iBandSize = *pBandCount;
        pSrcBand = new int[iBandSize];
        pDstBand = new int[iBandSize];

        for (int i=0; i<iBandSize; i++)
        {
            pSrcBand[i] = pBandIndex[i];
            pDstBand[i] = i+1;
        }
    }
    else
    {
        iBandSize = iBandCount;
        pSrcBand = new int[iBandSize];
        pDstBand = new int[iBandSize];

        for (int i=0; i<iBandSize; i++)
        {
            pSrcBand[i] = i+1;
            pDstBand[i] = i+1;
        }
    }
    
    void *hTransformArg = NULL, *hGenImgPrjArg = NULL;
    hTransformArg = hGenImgPrjArg = GDALCreateGenImgProjTransformer2((GDALDatasetH) pDSrc, (GDALDatasetH) pDDst, NULL);
    if (hTransformArg == NULL)
    {
        if(pProgress != NULL)
            pProgress->SetProgressTip("转?换?参?数?错?误?！?");

        GDALClose((GDALDatasetH) pDSrc);
        GDALClose((GDALDatasetH) pDDst);
        return RE_PARAMERROR;
    }
    
    GDALTransformerFunc pFnTransformer = GDALGenImgProjTransform;
    GDALWarpOptions *psWo = GDALCreateWarpOptions();

    psWo->papszWarpOptions = CSLDuplicate(NULL);
    psWo->eWorkingDataType = dataType;
    psWo->eResampleAlg = eResample;

    psWo->hSrcDS = (GDALDatasetH) pDSrc;
    psWo->hDstDS = (GDALDatasetH) pDDst;

    psWo->pfnTransformer = pFnTransformer;
    psWo->pTransformerArg = hTransformArg;

    psWo->pfnProgress = GDALProgress;
    psWo->pProgressArg = pProgress;

    psWo->nBandCount = iNewBandCount;
    psWo->panSrcBands = (int *) CPLMalloc(iNewBandCount*sizeof(int));
    psWo->panDstBands = (int *) CPLMalloc(iNewBandCount*sizeof(int));
    for (int i=0; i<iNewBandCount; i++)
    {
        psWo->panSrcBands[i] = pSrcBand[i];
        psWo->panDstBands[i] = pDstBand[i];
    }

    RELEASE(pSrcBand);
    RELEASE(pDstBand);

    GDALWarpOperation oWo;
    if (oWo.Initialize(psWo) != CE_None)
    {
        if(pProgress != NULL)
            pProgress->SetProgressTip("转?换?参?数?错?误?！?");

        GDALClose((GDALDatasetH) pDSrc);
        GDALClose((GDALDatasetH) pDDst);

        return RE_PARAMERROR;
    }

    oWo.ChunkAndWarpImage(0, 0, iNewWidth, iNewHeight);
    
    GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
    GDALDestroyWarpOptions( psWo );
    GDALClose((GDALDatasetH) pDSrc);
    GDALClose((GDALDatasetH) pDDst);

    if(pProgress != NULL)
        pProgress->SetProgressTip("重?采?样?完?成?！?");

    return RE_SUCCESS;
}

    PS：在使用Windows Live Writer来写博客，使用VSPaste插件粘贴代码的时候，发现会在汉字后面加一个“？”号，我懒得修改了，大家就凑合看看吧！

箜_Kong
关注
  • 4
    点赞
  • 踩
  • 37
    收藏
    觉得还不错? 一键收藏
  • 77
    评论
专栏目录
java结合gdal对大图片进行重采样
01-03
java结合gdal对大图片进行重采样，代码主要功能为制作缩略图，可以对较大的图片进行压缩处理，解决java不能处理太大文件的问题
gdal 投影转换 重采样方法
08-22
gadl 基本用法 重采样 重投影函数 封装成动态链接库的形式
基于gdal的遥感影像重采样
如果想成为中心，那么就到中心去吧。
11-20 1515
前言 遥感影像的重采样，不同于一般的图像，需要考虑空间参考。 代码 from osgeo import gdal, gdalconst import os def resampling(source_file, target_file, scale=5.): """ 影像重采样 :param source_file: 源文件 :param target_file: 输出影像 :param scale: 像元缩放比例 :return: """
【GDAL-Python】4-在Python中使用GDAL进行投影转换、影像重采样以及裁剪栅格数据
最新发布
weixin_44702962的博客
04-19 670
（1）在本教程中，将介绍如何在 Python 中使用 gdalwarp 将栅格数据重新投影到不同的坐标参考系，影像重采样（改变影像分辨率）、以及根据矢量 shapefile文件将影像对应部分进行裁剪； （2）技术路线：在python脚本中使用gdal.Warp，不同的功能对应不同的可选项参数1）坐标参考系改变对应：gdal.Warp(dstSRS = “EPSG:4326”)目标参考系，采用EPSG方式。
GDAL影像重采样
SUN 的博客
06-30 5998
学习过程中，用到重采样，感觉不错，转载过来供大家学习，原文地址：http://blog.csdn.net/giselite/article/details/7792620 #include "gdal_priv.h" #include "ogrsf_frmts.h" #include "gdalwarper.h" /***
【实用的GDAL】使用 gdalwarp 影像裁剪和重采样
weixin_43608276的博客
12-16 5835
使用gdalwarp对影像进行裁剪和重采样
python gdal 重采样_用GDAL对影像重采样的一些要点
weixin_39709674的博客
12-17 389
1、RasterIO只能使用最临近插值法做影像金字塔时，开始用的GDALDataSet的RasterIO通过控制行列宽度来抽层，代码写起来很简洁，奇怪一直没见设置重采样参数的地方，急于实现功能，没有怎么深究，基本功能全都做好后发现分层质量很差。文档里边没找到资料介绍这些，查了很久在Maillist上发现RasterIO() level resampling should all be "neare...
用GDAL对影像重采样的一些要点
weixin_33698043的博客
04-27 191
1、RasterIO只能使用最临近插值法 做影像金字塔时，开始用的GDALDataSet的RasterIO通过控制行列宽度来抽层，代码写起来很简洁，奇怪一直没见设置重采样参数的地方，急于实现功能，没有怎么深究，基本功能全都做好后发现分层质量很差。文档里边没找到资料介绍这些，查了很久在Maillist上发现RasterIO() level resampling should all be "ne...
【Python&RS】GDAL对栅格数据重采样/对齐栅格
Remote Sensing迷途小书童
05-26 2154
GDAL栅格数据重采样/对齐栅格
GDAL Warp重投影重采样
weixin_42069384的博客
11-04 1952
对于栅格数据的重投影可以使用Warp来实现 Warp(destNameOrDestDS, srcDSOrSrcDSTab,**kwargs) destNameOrDestDS是数据集变量或者新保存的目标文件地址。 srcDSOrSrcDSTab可以是源数据集或文件名，即被重投影的数据集，也可以是包含了若干数据集或文件名的一个array。 kwargs关键字参数是可以是一个WarpOptions对象，也可以是一个字符串数组，包含了WarpOptions的若干条件，具体内容与WarpOptions的参数一致。如
gdal应用-gdalwarp
弗里曼的小伙伴
03-31 6734
gdal应用-gdalwarp [email protected] 2015年1月16日 1 gdalwarp:变形工具。包括投影、拼接、及相关的变形功能。此工具功能强大，但效率不高，使用时注意 gdalwarp [--help-general] [--formats]     [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]   
Qt GDAL C++ 遥感图像显示 处理
05-05
软件包括了显示遥感图像，及鼠标拖动，随鼠标位置缩放，及部分栅格图像处理功能(图像锐化，平滑等)，和显示矢量图像的功能。
使用gdal来显示图像
01-04
使用gdal来显示图像，了解gdal如何使用，如何分块读取，显示
GDAL重采样与裁剪图像示例
豆芽菜
03-09 488
GDAL重采样，可以通过写文件时改变图像尺寸和geo_transformes的分辨率信息实现。GDAL裁切先计算索引，通过索引获取裁切后的图像，同时重新计算geo_transformes里的左上角点坐标。
基于python和GDAL实现遥感影像重采样（改变分辨率）
热门推荐
做GIS梦的人
05-23 1万+
很多时候都需要改变已经影像的分辨率，这里自己动手研究了一下相关原理，并进行了实现，以后可以很方便地改变影像的分辨率。 影像重采样核心 重采样的核心是影像的坐标范围不变，改变影像像元的大小，来实现像元个数的增减，即分辨率的改变。 像元的面积*像元个数=固定值 像元大小变为原来的1/2时，影像的像元数量变为原来的四倍。 代码实现 from osgeo import gdal, gdalconst import os import numpy as np def resampli..
Python地学分析 — GDAL对遥感影像重采样
一叶知秋的博客
12-12 8047
欢迎关注博主的微信公众号：“智能遥感”。 该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。 本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）： https://github.com/xbr2017 CSDN也在同步更新： https://blog.csdn.net/XBR_2014 “对遥感影像重采样是遥感...
GDAL源码剖析（十二）之GDAL Warp API使用说明
weixin_33896726的博客
05-25 728
一、简介 本文原文地址：http://www.gdal.org/warptut.html GDAL Warp API（在文件gdalwarper.h中定义）是一个高效的进行图像变换的接口。主要由几何变换函数（GDALTransformerFunc），多种图像重采样方式，掩码操作选项等组成。这个接口可以对很大的图像进行处理。 下面说明示例让你如何在程序中使用变换API。首先假定你已经熟...
java版GDAL学习之路—影像重采样
qq_41605733的博客
04-14 1402
目录整体思路代码实现 整体思路 对栅格影像进行重采样用到的主要方法是gdal的内置方法gdal.ReprojectImage()，采样后的数据需要新生成一个栅格数据集进行保存 思路： 载入原始栅格，读取相关信息 包括投影信息，仿射变换六参数，波段信息，数据类型，驱动名称，影像分辨率，影像宽、高等 新建输出栅格数据集 根据获取到的原始栅格的相关信息，以及重采样的因子（采样前分辨率/采样后分辨率），新建一个栅格数据集 进行重采样 gdal.ReprojectImage(dataset, outD
C++ gdal 重采样
07-28
您好！对于使用C++和GDAL库进行重采样的问题，您可以按照以下步骤进行操作： 1. 首先，确保您已经安装了GDAL库，并将其包含在您的C++项目中。您可以从GDAL官方网站（https://gdal.org/）下载并安装GDAL。 2. 在您的C++代码中，包含GDAL的头文件： ```cpp #include <gdal_priv.h> ``` 3. 创建一个GDAL数据集对象，并打开源图像文件： ```cpp GDALAllRegister(); // 注册所有的GDAL驱动 // 打开源图像文件 GDALDataset* srcDataset = (GDALDataset*)GDALOpen("path_to_source_image", GA_ReadOnly); if (srcDataset == nullptr) { // 处理打开源图像文件失败的情况 } ``` 4. 定义目标图像的宽度、高度和波段数目，并创建一个新的GDAL数据集对象： ```cpp int targetWidth = 800; // 目标图像宽度 int targetHeight = 600; // 目标图像高度 int targetBands = srcDataset->GetRasterCount(); // 目标图像波段数目与源图像相同 // 创建新的目标图像数据集 GDALDriver* memDriver = GetGDALDriverManager()->GetDriverByName("MEM"); GDALDataset* targetDataset = memDriver->Create("", targetWidth, targetHeight, targetBands, GDT_Byte, nullptr); if (targetDataset == nullptr) { // 处理创建目标图像数据集失败的情况 } ``` 5. 使用GDAL库提供的重采样方法，将源图像数据写入目标图像数据集中： ```cpp GDALResampleAlg resampleMethod = GRA_Bilinear; // 重采样方法，这里使用双线性插值 // 重采样 for (int band = 0; band < targetBands; ++band) { GDALRasterBand* srcBand = srcDataset->GetRasterBand(band + 1); // 获取源图像波段 GDALRasterBand* targetBand = targetDataset->GetRasterBand(band + 1); // 获取目标图像波段 // 执行重采样 GDALReprojectImage(srcBand, nullptr, targetBand, nullptr, resampleMethod); } ``` 6. 完成重采样后，您可以将目标图像保存到磁盘上： ```cpp GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff"); // 指定保存为GeoTIFF格式 driver->CreateCopy("path_to_target_image", targetDataset, 0, nullptr, nullptr, nullptr); ``` 7. 最后，别忘记释放使用的内存和关闭数据集： ```cpp GDALClose(srcDataset); // 关闭源图像数据集 GDALClose(targetDataset); // 关闭目标图像数据集 ``` 这些步骤提供了一个基本的框架，您可以根据您的具体需求进行修改和扩展。希望能对您有所帮助！如果您还有其他问题，请随时提问。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 77
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值